
Mitochondrial DNA Part A 

DNA Mapping, Sequencing, and Analysis 


ISSN: 2470-1394 (Print) 2470-1408 (Online) Journal homepage: http://www.tandfonline.com/loi/imdn21 


Taylor &. Francis 

Taylor & Francis Croup 


Phylogeography and spatial structure of the 
lowland tapir (Tapirus terrestris, Perissodactyla: 
Tapiridae) in South America 

Manuel Ruiz-Garcfa, Catalina Vasquez, Sergio Sandoval, Franz Kaston, Kelly 
Luengas-Villamil &Joseph MarkShostell 


To cite this article: Manuel Ruiz-Garcfa, Catalina Vasquez, Sergio Sandoval, Franz Kaston, 

Kelly Luengas-Villamil & Joseph Mark Shostell (2016) Phylogeography and spatial structure of the 
lowland tapir (Tapirus terrestris, Perissodactyla: Tapiridae) in South America, Mitochondrial DNA 
Part A, 27:4, 2334-2342, DOI: 10.3109/19401736.2015.1022766 

To link to this article: http://dx.d 0 i. 0 rg/l 0.3109/19401736.2015.1022766 


Published online: 22 May 2015. 


Submit your article to this journal GF 


Iilll Article 


views: 71 


5l View related articles GF 


View Crossmark data GF 


Full Terms & Conditions of access and use can be found at 
http://www.tandfonline.com/action/journallnformation?journalCode=imdn21 


r 


Download by: [Pontificia Universidad Javeria], [Manuel Ruiz-Garcia] 


Date: 19 July 2017, At: 10:51 



















DNA 


http://informahealthcare.com/mdn 

ISSN: 2470-1394 (print), 2470-1408 (electronic) 

Mitochondrial DNA Part A, 2016; 27(4): 2334-2342 
© 2015 Informa UK Ltd. DOI: 10.3109/19401736.2015.1022766 


informa 

healthcare 


SHORT COMMUNICATION 

Phylogeography and spatial structure of the lowland tapir 
[Tapirus terrestris, Perissodactyla: Tapiridae) in South America 
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Joseph Mark Shostell 4 
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and 4 Department of Biology, Penn State University-Fayette, Uniontown, Pennsylvania, USA 


Abstract 

We sequenced the mitochondrial cytochrome b gene of 141 lowland tapirs ( Tapirus terrestris) - 
representing the largest geographical distribution sample of this species studied across of 
South America to date. We compare our new data regard to two previous works on population 
structure and molecular systematics of T. terrestris. Our data agree with the Thoisy et al.'s work 
in (1) the Northern Western Amazon basin was the area with the highest gene diversity levels in 
T. terrestris, being probably the area of initial diversification; (2) there was no clear association 
between haplogroups and specific geographical areas; (3) there were clear population 
decreases during the last glacial maximum for the different haplogroups detected, followed by 
population expansions during the Holocene; and (4) our temporal splits among different 
T. terrestris haplogroups coincided with the first molecular clock approach carried out by these 
authors (fossil calibration). Nevertheless, our study disagreed regard to other aspects of the 
Thoisy et al.'s claims: (1) meanwhile, they detected four relevant clades in their data, we put 
forward six different relevant clades; (2) the Amazon River was not a strong barrier for 
haplotype dispersion in T. terrestris; and (3) we found reciprocal monophyly between T. terrestris 
and T. pinchaque. Additionally, we sequenced 42 individuals (T. terrestris, T. pinchaque, T. bairdii, 
and the alleged "new species", T. kabomani) for three concatenated mitochondrial genes ( Cyt-b, 
COI, and COII) agreeing quite well with the view of Voss et al., and against of the claims of 
Cozzuol et al. Tapirus kabomani should be not considered as a full species with the results 
obtained throughout the mitochondrial sequences. 
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Introduction 

In the Neotropics, only a few large mammalian species have been 
intensively studied in regard to their respective genetics and 
phylogeographic relationships (jaguar, Panthera onca ; Eizirik 
et al., 2001; Ruiz-Garcia, 2001; Ruiz-Garcia et al., 2003, 2006, 
2007, 2013; Andean bear, Tremarctos ornatus ; Ruiz-Garcia, 2003, 
2007, 2013; Ruiz-Garcia et al., 2003, 2005; the pink river dolphin, 
Inia geoffrensis and I. boliviensis , Banguera-Hinestroza et al., 
2002; Martmez-Aguero et al., 2006, 2010; Ruiz-Garcia, 2010a,b; 
Ruiz-Garcia et al., 2008). Similarly, the lowland tapir ( Tapirus 
terrestris ), the largest terrestrial herbivore in the South America, 
deserves to be studied more deeply of what have been studied so 
far from a molecular population genetics point of view. The aim 
of our study is to fill this gap in data. 

The Tapiridae family (Gray 1821), within the order 
Perissodactyla, is currently composed of one unique genus, 


Correspondence: Manuel Ruiz-Garcia, Laboratorio de Genetica de 
Poblaciones Molecular y Biologia Evolutiva, Departamento de 
Biologia, Facultad de Ciencias, Unidad de Genetica, Pontificia 
Universidad Javeriana, Bogota DC, Colombia. E-mail: mruizgar@ 
yahoo.es, mruiz@javeriana.edu.co 


Tapirus (Brisson 1762). Around 20 different Tapirus species are 
recognized for North America, Europa, and Asia in the fossil 
record (Hulbert, 2010). Currently, there are four species of 
Tapirus ; two exclusively present in South America (T. terrestris 
and T. pinchaque ), one from Southern Mexico, Central America, 
and the Pacific coasts of Colombia (T. bairdii) and another in Asia 
( T\ indicus). 

From a molecular genetics point of view, only a few works 
have been published with this genus. Ashley et al. (1996) and 
Norman & Ashley (2000) analyzed the genetic relationships 
among the Tapirus species by using mitochondrial Cytochrome 
Oxidase subunit II ( COII) and 12S rRNA gene sequences with a 
few individuals. More recently, the first molecular phylogeo¬ 
graphic study of T. terrestris analyzed 45 specimens at the 
mitochondrial cytochrome b gene {Cyt-b) (Thoisy et al., 2010). 
Ruiz-Garcia et al. (2012) showed that T. pinchaque did not show 
any significant spatial structure whereas T. bairdii showed a very 
significant spatial structure, due to isolation by distance. Cozzuol 
et al. (2013) claimed the existence of a new tapir species to 
compare Cyt-b gene sequences of four Amazon tapir individuals 
(two Brazilian animals and two Colombian specimens sampled by 
the first author of the present study) with regard to the 45 Cyt-b 
gene sequences published by Thoisy et al. (2010). Voss et al. 
(2014), however, put in doubt the validity of T. kabomani as a 
different species from T. terrestris. 
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We used sequences of the Cyt-b gene, which is commonly used 
among the molecular markers relevant for phylogeography, 
biosystematics, and genetic structure studies in mammal popula¬ 
tions (Patton et al., 2000). We sequenced 141 individuals of 
T. terrestris, including more geographical areas and individuals 
than any previous study (Cozzuol et al., 2013; Thoisy et al., 2010) 
to determine levels of gene diversity, spatial genetic structure, and 
demographic changes of T. terrestris in South America. 
Additionally, 42 tapir individuals from T. terrestris , T. pinchaque , 
T. bairdii , and the alleged new species, T. kabomani were 
sequenced at the Cyt-b , COI, and COII genes to compare with the 
tree showed by Cozzuol et al. (2013) and to show evidence in 
favor or against of the alleged “new species” T. kabomani. 

Methods 

Samples 

We sequenced a total of 141 individuals of T. terrestris at the 
Cyt-b gene from nine South American countries including 
Colombia (42 animals), Venezuela (six individuals), French 
Guiana (11 individuals), Ecuador (seven animals), Peru (30 
animals), Bolivia (11 animals), Brazil (24 animals), Paraguay 
(one animal), and Argentina (two animals) (Figure 1). 
Additionally, one animal from the Barcelona Zoo (Spain), five 
animals from the Cincinnati Zoo (Cincinnati, OH) and one animal 
of unknown origin were also analyzed. These last seven animals 
were not included in those analyses where the origin of the 
animals was necessary. Thirty T. pinchaque (12 haplotypes), 
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30 T. bairdii (six haplotypes), and one T. indicus samples obtained 
by the first author were also employed. For the question of 
T. kabomani , 42 individuals were analyzed: 17 T. terrestris (one 
from Argentina, two from Bolivia, three from Brazil, five from 
Ecuador, three from Peru and three from the Cincinnati Zoo), five 
‘T. kabomani ” (the two Brazilian individuals reported by 
Cozzuol et al., 2013, and three “classical morphological” 
T. terrestris sampled by us but with haplotypes related to 
T. kabomani ; one from San Martin de Amacayacu, Amazon River, 
Colombia, one from the Mazan River, a tributary from the Napo 
River, in the Peruvian Amazon and one from near to Tena, upper 
Napo River, at the Ecuadorian Amazon), 12 T. pinchaque (three 
Colombian and nine Ecuadorian individuals), and nine T. bairdii 
(all them from the Darien area in Panama). All these 42 
individuals were sequenced for Cyt-b (1140 base pairs, bp), 
COI (640 bp), and COII (726 bp) genes summing up 2506 bp. 

Molecular analyses 

DNA from teeth, bones, muscle, and skins were obtained with the 
phenol-chloroform procedure (Sambrock et al., 1989), whereas 
DNA samples from hair and blood were obtained with 10% 
Chelex® 100 resin (Walsh et al., 1991). Amplifications for Cyt-b 
gene were achieved using primers designed for perissodactyls 
(Tougard et al., 2001). The amplification conditions employed 
followed Ruiz-Garcia et al. (2012). For the COI and COII genes, 
the amplification conditions followed Ashley et al. (1996) and 
Hebert et al. (2003, 2004). All amplifications, including positive 



N. Locality 
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Figure 1. Map of the geographical localities where 141 lowland tapirs were sampled to sequence the Cyt-b gene. 
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and negative controls, were checked in 2% agarose gels, 
employing the molecular weight marker 0X174 DNA digested 
with Hind III, Hinf I genes, and HyperLadder IV. Those samples 
that amplified were purified using membrane-binding spin 
columns (Qiagen, Valencia, CA). The double-stranded DNA 
was directly sequenced in a 377A (ABI, Vernon, CA) automated 
DNA sequencer. The samples were sequenced in both directions 
and all the samples were repeated to ensure sequence accuracy. 
The sequences were deposited in Genbank (accession nos. GQ 
259910-1 and GQ 259954-1). 

Data analyses 

Genetic diversity and heterogeneity analyses 

The sequences were edited and aligned with BioEdit Sequence 
Alignment Editor (Hall, 2004) and DNA Alignment (Fluxus 
Technology Ltd, Charlotte, NC). 

The statistics employed to determine the genetic diversity at 
the Cyt-b gene for T. terrestris were the number of polymorphic 
sites (S), the number of haplotypes (H), the haplotypic diversity 
(H d ), the nucleotide diversity (7r), the average number of 
nucleotide differences ( k ), and the 0 statistic by sequence. 
Different tests were carried out to measure genetic heterogeneity, 
and possible gene flow estimates, among the different T. terrestris 
groups, were considered. These tests were those of Hudson et al. 
(1992) (H s t, K s T , X ST *, Z, and Z*), Snn test, and the Chi-square 
test on the haplotypic frequencies with permutation tests using 
10,000 replicates. We also estimated the G ST statistic from the 
haplotypic frequencies and the y S T, N ST , and F ST statistics 
(Hudson et al., 1992) from the nucleotide sequences. All these 
statistics were obtained by using the software DNAsp 5.10 
(Florida Institute of Technology, Vero Beach, FL) (Librado & 
Rozas, 2009). 

Phylogenetic analyses 

To determine phylogenetic relationships and temporal splits 
among the different haplotype clades within T. terrestris at the 
Cyt-b gene, a Bayesian tree was performed using a GTR (General 
Time Reversible) model of nucleotide substitution. We used the 
BEAST v. 1.6.2 program (Viacom International Inc., New York, 
NY) for this analysis, assuming a Yule speciation model and a 
relaxed molecular clock with an uncorrelated log-normal rate of 
distribution (Drummond et al., 2006). Results from the two 
independent runs (40,000,000 generations with the first 4,000,000 
discarded as burn-in) were combined with LogCombiner vl.6.2 
software (Gene Code Corporation, Ann. Arbor, MI) (Viacom 
International Inc., New York, NY) (Drummond & Rambaut, 
2007). Posterior probability values of each node were estimated, 
which provided an assessment of the degree of support of each 
node on the tree. We employed a calibration point of 18 ±1 
millions of year ago (MYA) between the ancestor of T. indicus 
and the ancestor of the Neotropical tapirs and the calibration point 
between the ancestor of T. bairdii and the ancestor of the two 
South American species, T. terrestris + T. pinchaque was taken as 
9.5 ±0.5 MYA following the paleontological and molecular 
results of different authors (Ashley et al., 1996; Colbert, 2005; 
Hulbert, 1995, 2010; Hulbert & Wallace, 2005; Hulbert et al., 
2009; Norman & Ashley, 2000; Ruiz-Garcia et al., 2012). No 
temporal split priors, or monophylia, between T. terrestris and 
T. pinchaque , were imposed by the authors, being the temporal 
split between these taxa estimated by the software. 

To obtain proofs in favor or against of T. kabomani , we 
obtained a maximum likelihood tree (Felsenstein, 1981) for the 
three mitochondrial genes concatenated to compare with the same 
tree from Cozzuol et al. (2013). 


Demographic changes 

A Bayesian skyline plot (BSP) was obtained by means of the 
BEAST v. 1.6.2 and Tracer vl.5 software (Viacom International 
Inc., New York, NY). The Coalescent-Bayesian skyline option in 
the tree priors was selected with five steps and a piecewise- 
constant skyline model with 40,000,000 generations (the first 
4,000,000 discarded as burn-in). A stepwise (constant) Bayesian 
skyline variant was selected with the maximum time as the upper 
95% higher posterior density (HPD). 

Spatial genetic analysis on T. terrestris 

A spatial autocorrelation analysis (Sokal & Oden, 1978a,b; Sokal 
& Wartenberg, 1983) was applied. Coordinates were given to each 
one of the individuals sequenced and the AIDA procedure 
implemented by Bertorelle and Barbujani (1995) was applied. We 
calculated the autocorrelation coefficient II and their respective 
correlograms (equivalent to the traditional Moran’s I index; 
Moran, 1950). To connect the individuals sequenced within each 
distance class, the Gabriel-Sokal network (Gabriel & Sokal, 1969; 
Matula & Sokal, 1980) was employed. 

Results and discussion 

The 141 individuals of T. terrestris we analyzed showed high 
levels of genetic diversity. These gene diversity statistics were 
also calculated in three different geographical levels: within the 
haplogroups found, by geographical regions and by being north or 
south of the Amazon River (Table 1). All the nucleotide 
diversities are similar within each haplogroup, with both the 
North and the Amazon III haplogroups, the lowest. By geograph¬ 
ical region, the Northern Western Amazon yielded the highest 
nucleotide diversity, whereas the French Guyana and the Southern 
South America region were those which presented the lowest 
nucleotide diversities. If we analyze all the individuals sampled to 
the north and south of the Amazon River, the north set showed a 
higher nucleotide diversity. 

Our results showed similar levels of gene diversity to those 
estimated by Thoisy et al. (2010) (H d — 0.988 and 
7r — 0.0093 ± 0.0004), although in our sample (141 versus 45 
individuals), we detected a considerable higher number of 
haplotypes (80 versus 35). The highest levels of gene diversity 
were located in the lowland tapir population from the Western 
northern Amazon basin. This agrees quite well with the fact that 
the most ancestral and original T. terrestris population was this 
one. This result ratifies findings by Thoisy et al. (2010).The 
genetic heterogeneity among different tapir sets can be seen at 
Table 2. The classification of haplogroups is based on the highest 
levels of genetic heterogeneity among the different tapir sets. 

The seven genetic heterogeneity tests indicated significant 
findings. Additionally, the F ST statistic was considerably large 
(F ST = 0.716). The seven heterogeneity tests were significantly 
different by geographical regions, but the values of F ST were 
considerably smaller than in the previous case (F ST = 0.279). This 
is a consequence of different haplotypes belonging to different 
haplogroups and being intermixed within the same geographical 
area. Finally, we compared the genetic heterogeneity between the 
animals sampled north and south of the Amazon River. Although 
the seven genetic heterogeneity tests were significant, the value of 
F ST (= 0.045) was small, which could support that the Amazon 
River has not been a complete geographical barrier to the 
dispersion of the tapir. 

Thoisy et al. (2010) affirmed that the Amazon River acted as a 
barrier to gene flow. However, our results showed that the 
Amazon River was only a partial barrier for haplotype dispersion 
in T. terrestris (Nm =11-33). 


DOi: 10.3109/19401736.2015.1022766 Phylogeography of the South American lowland tapir 2337 


Table 1. Gene diversity statistics for Tapirus terrestris by mitochondrial lineages, geographic regions, and different Amazon River banks at the 
Cyt-b gene. 



S 

NH 

H d 

7 r 

K 

0 per sequence 

Total sample studied of Tapirus terrestris 

107 

80 

0.984 ±0.003 

0.0114 ±0.0003 

10.335 ±4.743 

19.376 ±4.739 

Mitochondrial lineages 

Amazon I 

14 

10 

0.949 ±0.051 

0.0038 ±0.0019 

3.436 ±1.876 

4.51 ±2.024 

Amazon II 

40 

20 

0.937 ±0.020 

0.0039 ±0.0005 

3.563 ±1.849 

6.583 ±2.232 

Amazon III 

30 

20 

0.926 ±0.025 

0.0043 ±0.0006 

3.937 ±2.017 

7.186 ±2.442 

North 

27 

18 

0.921 ±0.034 

0.0051 ±0.0005 

4.065 ±2.082 

6.704 ±2.362 

South 

20 

10 

0.867 ±0.079 

0.0040 ±0.0008 

3.642 ±1.946 

6.027 ±2.479 

Geographical regions 

Northern Colombia and Venezuela 

27 

14 

0.948 ±0.031 

0.0065 ±0.0002 

5.943 ±2.953 

7.505 ±2.835 

Guyana region 

10 

4 

0.600 ±0.154 

0.0033 ±0.0011 

2.982 ±1.685 

3.414 ±1.664 

Upper Northern Amazon 

69 

40 

0.973 ±0.009 

0.0129 ±0.0001 

11.779 ±5.402 

14.545 ±4.140 

Southern Amazon (Southern Peru and Bolivia) 

33 

15 

0.910 ±0.044 

0.0091 ±0.0002 

8.263 ±3.962 

8.740± 3.131 

Southern South America 

11 

5 

0.806 ±0.120 

0.0034 ±0.0002 

3.779 ±1.779 

4.545 ±2.010 

Amazon River 

Northern Amazon river bank 

86 

58 

0.980 ±0.005 

0.0121 ±0.0002 

10.969 ±5.029 

16.643 ±4.361 

Southern Amazon river bank 

51 

24 

0.945 ±0.022 

0.0089 ±0.0004 

8.154 ±3.863 

11.990 ±3.794 


The statistics estimated were the number of polymorphic sites (S), the number of haplotypes (NH), the haplotypic diversity (H d ), the nucleotide 
diversity ( 7 r), the average number of nucleotide differences (AT), and the 0 statistic by sequence. 


Table 2. Genetic heterogeneity statistics for Tapirus terrestris by 
mitochondrial lineages, geographic regions, and by the Amazon River 
banks at the Cyt-b gene. 


Genetic differentiation estimated 

P 

Gene flow 

Mitochondrial lineages 




X 2 = 696.063, df= 390 

0.0000* 

G st = 0.071 

Nm = 6.54 

H st = 0.0609 

0.0000* 

7st = 0.6384* 

Nm = 0.28 

^st = 0.6264 

0.0000* 

Vst = 0.7177* 

Nm = 0.20 

^st* = 0.3752 

0.0000* 

F st = 0.7160* 

Nm = 0.20 

Z s = 1403.773 

0.0000* 



Z s * = 6.8447 

0.0000* 



5^ = 0.9857 

0.0000* 



Geographical regions 




X 2 = 457.206, df= 280 

0.0000* 

G st = 0.0732 

Nm = 6.33 

H st = 0.0641 

0.0000* 

7st = 0.1509* 

Nm = 2.81 

^st = 0.1333 

0.0000* 

N st = 0.2789* 

Nm=1.29 

^ ST * = 0.1215 

0.0000* 

F st = 0.2794* 

Nm=1.29 

Z s = 3896.878 

0.0000* 



Z s * = 7.6261 

0.0000* 



Sim = 0.7661 

0.0000* 



Amazon River 




X 2 = 120.256 df == 77 

0.0000* 

G st = 0.0148 

Nm = 33.30 

H st = 0.0131 

0.0000* 

7st = 0.0243 

Nm = 20.09 

^ ST = 0.0180 

0.0000* 

N st = 0.0449 

Nm= 10.62 

^st* = 0.0186 

0.0000* 

F st = 0.0450 

Nm= 10.60 

Z s = 4727.594 

0.0000* 



Z s * = 8.1055 

0.0000* 



5^ = 0.8377 

0.0000* 




*Significant probability (p< 0.0001). 


Six well-differentiated mitochondrial lineages in T. terrestris , 
with significant genetic heterogeneity among them, support the 
fact that different events occurred, which fragmented and isolated 
these lineages, in the past. However, after these isolation events, 
these different mitochondrial lineages again migrated, expanded, 
and mixed in sympatrical areas. It is for this reason that there are 
several mitochondrial lineages in the same geographical area 
today. For instance, six out of six mitochondrial haplogroups were 
found in the Colombian Amazon and four out of six were found in 
the Northern Peruvian Amazon. 

We detected two new Amazon haplogroups which go 
unnoticed by Thoisy et al. (2010): the Amazon 0 (clade very 
relevant because coincides with that Cozzuol et al., 2013 claimed 


as a different species, T. kabomani) and another Amazonian 
haplogroup (Amazon III). It is very easy to determine the 
correspondence of our remainder haplogroups with the four clades 
of Thoisy et al. (2010) because we provide the major part of the 
samples analyzed in Thoisy et al. (2010), and some of them were 
also enclosed in this work: Amazon I is the clade II, Amazon II is 
the clade I, North is the clade III, and South is the clade IV. 

Phylogeny and phylogeography of Tapirus terrestris 

The Bayesian tree (Figure 2) showed six haplogroups within T. 
terrestris at the Cyt-b gene. The posterior probabilities (p) of 
these main haplogroups were elevated in the majority of the cases: 
Amazon 0 (p = 1), Amazon I (p = 1), Amazon II (p = 1), Amazon 
III (p = 1), North ( p = 0.84), and South (p = 0.23). Amazon 0 and 
Amazon I were the first two clades in which ancestors first 
diverged inside T. terrestris. Amazon 0 was composed by two 
exemplars sampled at the Colombian Amazon (San Martin de 
Amacayacu) and at the Mazan River (a tributary of the Napo 
River) in the northern Peruvian Amazon. These two individuals 
showed haplotypes related to T. kabomani sensu Cozzuol et al. 
(2013). Amazon I was composed of 13 individuals (Colombian 
Amazon, Northwestern Brazilian Amazon, Northern Peruvian 
Amazon and Ecuadorian Amazon). Amazon II was composed of 
40 individuals (three Colombian Amazon departments, 
Colombian Eastern Llanos, Northern Colombia, Northern and 
Southern Peruvian Amazon, Central Brazilian Amazon, 
Northwestern Brazilian Amazon, Bolivian Amazon and 
Ecuadorian Amazon). Amazon III was composed of 37 individ¬ 
uals (Colombian Amazon, Eastern Colombian Llanos, 
Northwestern Brazilian Amazon, Central Brazilian Amazon, 
Northern and Southern Peruvian Amazon, Bolivian Amazon, 
and French Guyana). South was composed of 16 individuals 
(Southern Brazil, Paraguay, Bolivia, and Central Brazilian 
Amazon). The last haplogroup was named North and it was 
composed of 32 individuals. It was constituted by individuals 
from (French Guyana, Eastern Colombian Llanos, Northern 
Colombia and Venezuela and some individuals from Argentina, 
Southern Peruvian Amazon, Colombian Amazon, and 
Northwestern Brazilian Amazon). 

We estimated the temporal splits of the main nodes in the 
Bayesian tree. The temporal split between the ancestors of 
T. indicus and the Neotropical tapirs was around 18.29 MYA. 
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I T. indicus 
T. bairdii 


T. terrestris 
Amazon II 


T. terrestris 
North 


T. terrestris 
South 


T. terrestris 
Amazon III 


T. terrestris 
Amazon 0 

T. terrestris 
Amazon I 


T. pinchaque 


Figure 2. Bayesian tree (Beast program) with the haplotypes obtained from 141 Tapirus terrestris individuals sequenced at the Cyt-b gene. First 
numbers in the nodes are the posterior probabilities and second numbers are divergence time splits in millions of years ago. This Bayesian tree also 
encloses one haplotype (one individual) of T. indicus from Thailand, 12 haplotypes (30 individuals) of T. pinchaque from Colombia and Ecuador, and 
six haplotypes (30 individuals) of T. bairdii from Mexico, Guatemala, Costa Rica, and Panama. 


The ancestor of T. bairdii diverged from the ancestor of 
T. pinchaque + T. terrestris around 9.6 MYA. The split between 
T. pinchaque and T. terrestris was estimated around 3.33 MYA 
with a haplotype diversification process within the first taxon 


around 2.54 MYA. The ancestors of the Amazon 0 + Amazon I 
haplogroups diverged from the ancestors of the other T. terrestris 
haplogroups around 2.77 MYA, whereas the diversifications 
within the other haplogroups occurred slightly more recently and 
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were similar to each other (Amazon III and South, 2.43 MYA; 
North, 2.18 MYA and Amazon II, 2.09 MYA). There are three 
periods of haplotype diversification within the Amazon I, II, and 
III haplogroups as well as in the North and South haplogroups. 
These three temporal splits occurred around 1.2-1.7 MYA, 
0.6-0.9 MYA, and 0.25-0.14 MYA, respectively. 

Two main results from a phylogenetic point of view were 
obtained by Thoisy et al. (2010). The first one is that their 
Bayesian analysis favored T. pinchaque paraphyly in relation to 
lowland tapir over reciprocal monophyly. Our tree did not agree 
with the Thoisy et al. (2010)’ conclusion. Our Bayesian tree 
clearly showed that there was reciprocal monophyly between 
T. terrestris and T. pinchaque. No constriction was imposed in 
this tree to force this reciprocal monophyly. Second, Thoisy et al. 
(2010) employed two approaches to estimate temporal splits 
within the clades of T. terrestris. In the first approach, they used 
fossil data to calibrate the molecular clock with the timing 
of the Rhinocerontidae and Tapiridae split (46.7 MYA) and 
another split within Rhinocerontidae (17.1 MYA). In the sec¬ 
ond approach, they employed a direct mutation rate at the Cyt-b 
gene for Perissodactyls and they obtained a very recent mito¬ 
chondrial diversification within T. terrestris in the mid to late 
Pleistocene (0.33 MYA). We employed a different approach to 
calibrate the molecular clock by using simultaneously mixed 
fossil and molecular data. Our split results agree considerably 
better with the first approach of Thoisy et al. (2010) than with the 
second one. 

Demographic changes 

Our analyses of BSP findings indicate several important points. 
For the total T. terrestris sample, there was a population 
contraction around 2.2-2.0 MYA. However, 1.9 MYA, there 
was a strong population expansion, which reached its “plateau” 
around 1.5 MYA. Since then, the overall tapir population has 
remained constant until today. By haplogroups, Amazon I began 
decreasing 5500 YA and continued to do so until around 1400 
YA. Since then, this population has expanded up until the last 500 
years when it was constant. Amazon II and III experienced steady 
decreasing population trends beginning 75,000 YA and ending 
12,000 YA in the case of Amazon II and 14,000 YA in the case of 
Amazon III. Both haplogroups have since increased up until the 
last 500 years when their size became constant. North experienced 
a constant decreasing population trend that started around 
14,000 YA and ended 3000 YA, after which the lineage population 
increased until today. South experienced a slight but constant 
decreasing population trend from 6000 YA until 2500 YA. This 
lineage has increased since that time until today. 

Amazon I suffered a population declination 5500 YA coincid¬ 
ing with one of the dry periods in the Holocene after the Optimum 
Climaticum (OP). This dry period was detected in the Amazon, 
Caqueta and lower Magdalena River basins as well as in Andean 
lagoons in Colombia and Peru (Thompson et al., 1995; Van der 
Hammen, 2001; Van der Hammen & Cleef, 1992). Amazon II and 
III revealed a population decrease 75,000 YA coinciding with the 
ending of the Eemiense inter-glacial period and the beginning of 
the upper Pleniglacial period (Van der Hammen, 1992). These 
two Amazonian haplogroups began to again expand around 
12,000 and 14,000 YA after LGM. In contrast, North decreased 
14,000 YA, agreeing with the massive extinction of mammals 
across the Earth including South America during LGM. This 
corresponds with the Younger Dry as (Dry as III), typical of 
Northern Europe and Scandinavia (Clapperton, 1993). Finally, 
South suffered a population declination around 6000 YA coincid¬ 
ing with the drier period between 7000 and 5500 YA that we 
commented on above. 
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Spatial genetic structure 

The AIDA analyses generated different correlograms. This 
analysis, applied to the total tapir population, indicated that 
the individuals separated by 300 km or less (first distance class, 
1 DC) were significantly more related than expected by random, 
whereas the animals separated by 700 and 1300 km showed a 
significant negative spatial autocorrelation (3 and 4 DC). After 
1300 km, the genetic relationships among the individuals were 
more randomly distributed following a “crazy quilt” pattern 
following the suggestion of Sokal & Oden (1978a). The 7 DC 
(1650-2150 km) was positively significant (which suggests 
genetic flow among animals separated by this distance interval), 
whereas between 2150 and 2600 km, the autocorrelation was 
significantly negative (which suggests some obstacle or some 
geographic barrier at this distance interval). Thus, this AIDA 
analysis showed genetics patches with a diameter of around 
1300 km. The AIDA analysis - taking all the individuals 
sampled in Colombia and Venezuela, showed a correlogram 
typically representing an isolation by distance pattern. The 
individuals within the 1 DC (less than 350 km) showed signifi¬ 
cant genetic resemblance. However, animals separated by 900- 
1800 km (3 and 4 DC) were significantly different. The AIDA 
analysis applied to the Western Amazon showed a unique 
distance class that was significantly positive (4 DC, 
490-690 km). Therefore, migration events seem to have occurred 
at these distances. This particular area of the Amazon seems to 
have less spatial genetic structure compared to other geograph¬ 
ical areas analyzed. Finally, an AIDA analysis was applied to 
those tapirs sampled at the southern Amazon. The correlogram 
showed two significant distance classes. The first was 1 DC 
(0-150 km), where the autocorrelation was positive. For the 
second one, 4 DC (780-900 km), the autocorrelation was nega¬ 
tive. Thus in this Southern Amazonian area, the presence of 
a monotonic cline was detected. 

There was no clear association between haplogroups and 
specific geographical areas - revealing a complex system of 
migration and colonization by the lowland tapir, which coincides 
with the Thoisy et al.’s (2010) claim to reject the null hypothesis 
of allopatric divergence. However, we provide much more robust 
evidence of this complex system of migration and colonization 
than the Thoisy et al.’s (2010) work throughout the autocorrel¬ 
ation analysis. The application of AIDA revealed the existence of 
genetics patches of a diameter around 1300 km for the overall 
tapir population, but with repetitive events of gene flow at long 
distances (up until 5000 km). These genetics patches agree quite 
well with the existence of the Pleistocene Refugia (Haffer, 1969, 
1997, 2008). The AIDA analysis applied to the Upper Western 
Amazon basin indicated the smallest degree of spatial structure. 
This could be related with the fact that this was the population 
with less evidence of population expansion, and with the highest 
gene diversity levels because this was probably the original and 
more demographically stable lowland tapir population in South 
America. This agree quite well with the Napo Refugium (or Napo 
island) both claimed by the Pleistocene Refugia hypothesis 
(Whithmore & Prance, 1987) and the Recent Amazon lake 
hypothesis (Nores, 1999, 2004). 

The question of “T. kabomani" 

Figure 3 shows the maximum likelihood tree for the three 
mitochondrial genes concatenated. There is no evidence that 
T. kabomani should be a full species, being only a haplogroup 
(clade) within T. terrestris. Clearly, to augment the number of 
T. pinchaque individuals sequenced with reference to the work of 
Cozzuol et al. (2013), T. pinchaque is more differentiated from 
T. terrestris than T. kabomani is. 
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Figure 3. Maximum likelihood tree with three concatenated mitochondrial genes ( Cyt-b + COI + COII) for 42 tapirs, including T. terrestris, 
T. kabomani, T. pinchaque, and T. bairdii. Tapirus kabomani is a clade within T. terrestris. There was reciprocal monophyly between T. terrestris and 
T. pinchaque. 


The question of the alleged new species, T. kabomani , is 
relevant for the knowledge of the phylogeography of T. terrestris 
because we provide proofs that T. kabomani is a haplogroup 
within T. terrestris. Cozzuol et al. (2013) affirmed the existence of 
T. kabomani because four sequences of two Colombian and two 
Brazilian tapirs were outside of the T. terrestris + T. pinchaque 
clade. However, our maximum likelihood tree showed that 


T. kabomani is more related to T. terrestris than T. pinchaque is. 
The inaccurate result obtained by Cozzuol et al. (2013) was 
probably related with the fact that these authors only analyzed five 
samples of T. pinchaque at the Cyt-b gene and only one sample at 
the Cyt-b + COI + COII genes. The very small T. pinchaque 
sample employed by Cozzuol et al. (2013) did not probably 
represent all the mitochondrial gene diversity of T. pinchaque and, 
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as the genetic differences between both T. terrestris and T. 
pinchaque are very small, by chance the no representative 
mitochondrial gene diversity of the T. pinchaque of the small 
sample of Cozzuol et al. (2013) was nested within the gene 
diversity of T. terrestris. Nevertheless, at the moment that we 
enlarged the T. pinchaque sample, this phenomenon disappeared. 
Henceforth, our tree, which contained the highest number of 
individuals of T. terrestris , T. kabomani , and T. pinchaque 
analyzed until date for the three concatenated mitochondrial genes 
employed by Cozzuol et al. (2013), agrees extremely well with the 
point of view of Voss et al. (2014) and against Cozzuol et al.’s 
(2013) claims. 

We believe more prudent to consider T. kabomani as a 
T. terrestris haplogroup until new data, employing nuclear and 
MHC gene sequences, have been generated for this taxon. 
Especially important in our view is to obtain karyotypes of this 
taxon to determine, or no, possible stasipatric speciation (due to 
the inexistence of geographical barriers between T. terrestris and 
the alleged T. kabomani ; White, 1968, 1978). 
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